### packages
library(cjoint)
library(foreign)
library(dplyr)
library(cregg)
library(ggplot2)
library(patchwork)
library(labelled)

########################
# Load and set up data #
########################
## France
# load  data
dat.e1fr <- read.csv("France - cjoint 1.csv")
# organize factor variable order
dat.e1fr$ideology <- factor(dat.e1fr$ideology, levels = c(0,1,2,3,4))
dat.e1fr$leadexp <- factor(dat.e1fr$leadexp, levels = c("Media","Business","Politics"))
dat.e1fr$selection <- factor(dat.e1fr$selection, levels = c("Leadership","Delegates","Members","All voters"))
dat.e1fr$status <- factor(dat.e1fr$status, levels = c("New party","Opposition","Government"))
dat.e1fr$women <- factor(dat.e1fr$women, levels = c("Below average","Above average"))

## Israel
## load recoded data 
dat.e1il <- read.csv("Israel - cjoint 1.csv")
## organize factor variable order
dat.e1il$ideology <- factor(dat.e1il$ideology, levels = c(0,1,2,3,4))
dat.e1il$leadexp <- factor(dat.e1il$leadexp, levels = c("Media","Business","Politics"))
dat.e1il$selection <- factor(dat.e1il$selection, levels = c("Leadership","Delegates","Members","All voters"))
dat.e1il$status <- factor(dat.e1il$status, levels = c("New party","Opposition","Government"))
dat.e1il$women <- factor(dat.e1il$women, levels = c("Below average","Above average"))

## Sweden
## load recoded data 
dat.e1sw <- read.csv("Sweden - cjoint 1.csv")
## organize factor variable order
dat.e1sw$ideology <- factor(dat.e1sw$ideology, levels = c(0,1,2,3,4))
dat.e1sw$leadexp <- factor(dat.e1sw$leadexp, levels = c("Media","Business","Politics"))
dat.e1sw$selection <- factor(dat.e1sw$selection, levels = c("Leadership","Delegates","Members","All voters"))
dat.e1sw$status <- factor(dat.e1sw$status, levels = c("New party","Opposition","Government"))
dat.e1sw$women <- factor(dat.e1sw$women, levels = c("Below average","Above average"))


#######################################
# estimation - preferences (Figure 3) #
#######################################
## estimate AMCE
out.amce.fr1 <- cj(data = dat.e1fr,
                   selected ~ ideology + selection + leadexp + status + women,
                   id = ~Response.ID,
                   estimate = "amce",
)
out.amce.il1 <- cj(data = dat.e1il,
                   selected ~ ideology + selection + leadexp + status + women,
                   id = ~Response.ID,
                   estimate = "amce",
)
out.amce.sw1 <- cj(data = dat.e1sw,
                   selected ~ ideology + selection + leadexp + status + women,
                   id = ~Response.ID,
                   estimate = "amce",
)
## plot outcomes
# define attribute order for graphs
out.amce.fr1$feature <- factor(out.amce.fr1$feature, 
                               levels = c("selection","ideology","leadexp","status","women"))
out.amce.il1$feature <- factor(out.amce.il1$feature, 
                               levels = c("selection","ideology","leadexp","status","women"))
out.amce.sw1$feature <- factor(out.amce.sw1$feature, 
                               levels = c("selection","ideology","leadexp","status","women"))
# define labels for graphs
attributes <- c("Ideological \nDistance","Candidae \nSelection","Leader \nExperience","Party \nStatus",
                "Female \nRepresentation") 
names(attributes) <- c("ideology","selection","leadexp","status","women")
# plot graphs
plot.amce.fr1 <- plot(out.amce.fr1, size=1.5, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
                      labeller = labeller(.rows = attributes)) +
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower", xmax = "upper"),  
                          size = 0.6, height = 0, na.rm = TRUE,
                          position = ggstance::position_dodgev(height = 0.75)) +
  ggplot2::scale_color_manual(values = rep("black", 17)) +
  ggplot2::xlim(-0.6,0.3) +
  ggplot2::theme(legend.position = "none") +
  ggplot2::ggtitle("France")
plot.amce.il1 <- plot(out.amce.il1, size=1.5, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
                      labeller = labeller(.rows = attributes)) +
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower", xmax = "upper"),  
                          size = 0.6, height = 0, na.rm = TRUE,
                          position = ggstance::position_dodgev(height = 0.75)) +
  ggplot2::scale_color_manual(values = rep("black", 17)) +
  ggplot2::xlim(-0.6,0.3) +
  ggplot2::theme(legend.position = "none", axis.text.y=element_blank()) +
  ggplot2::ggtitle("Israel")
plot.amce.sw1 <- plot(out.amce.sw1, size=1.5, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
                      labeller = labeller(.rows = attributes)) +
  ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower", xmax = "upper"),  
                          size = 0.6, height = 0, na.rm = TRUE,
                          position = ggstance::position_dodgev(height = 0.75)) +
  ggplot2::scale_color_manual(values = rep("black", 17)) +
  ggplot2::xlim(-0.6,0.3) +
  ggplot2::theme(legend.position = "none", axis.text.y=element_blank()) +
  ggplot2::ggtitle("Sweden")
plot.acme.all1 <- plot.amce.fr1 + plot.amce.il1 + plot.amce.sw1 
# save file (Figure 3)
ggsave(file = "Figure3.pdf", plot=plot.acme.all1, width=8, height=6.5)

######################################
# estimation - importance (Figure 4) #
######################################
### recode data - subset choices with cost in a given attribute compared to paired alternative
## lower ideological proximity
# france
dat.e1fr$ideology.int <- as.numeric(dat.e1fr$ideology)
dat.e1fr.lessideology <- dat.e1fr %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & ideology.int>=1 & lag(ideology.int)==0) |
      (task==dplyr::lead(task) & ideology.int>=1 & dplyr::lead(ideology.int)==0) |
      (task==lag(task) & ideology.int>=2 & lag(ideology.int)==1) |
      (task==dplyr::lead(task) & ideology.int>=2 & dplyr::lead(ideology.int)==1) |
      (task==lag(task) & ideology.int>=3 & lag(ideology.int)==2) |
      (task==dplyr::lead(task) & ideology.int>=3 & dplyr::lead(ideology.int)==2) |
      (task==lag(task) & ideology.int>=4 & lag(ideology.int)==3) |
      (task==dplyr::lead(task) & ideology.int>=4 & dplyr::lead(ideology.int)==3)
  )
# Israel
dat.e1il$ideology.int <- as.numeric(dat.e1il$ideology)
dat.e1il.lessideology <- dat.e1il %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & ideology.int>=1 & lag(ideology.int)==0) |
      (task==dplyr::lead(task) & ideology.int>=1 & dplyr::lead(ideology.int)==0) |
      (task==lag(task) & ideology.int>=2 & lag(ideology.int)==1) |
      (task==dplyr::lead(task) & ideology.int>=2 & dplyr::lead(ideology.int)==1) |
      (task==lag(task) & ideology.int>=3 & lag(ideology.int)==2) |
      (task==dplyr::lead(task) & ideology.int>=3 & dplyr::lead(ideology.int)==2) |
      (task==lag(task) & ideology.int>=4 & lag(ideology.int)==3) |
      (task==dplyr::lead(task) & ideology.int>=4 & dplyr::lead(ideology.int)==3)
  )
# Sweden
dat.e1sw$ideology.int <- as.numeric(dat.e1sw$ideology)
dat.e1sw.lessideology <- dat.e1sw %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & ideology.int>=1 & lag(ideology.int)==0) |
      (task==dplyr::lead(task) & ideology.int>=1 & dplyr::lead(ideology.int)==0) |
      (task==lag(task) & ideology.int>=2 & lag(ideology.int)==1) |
      (task==dplyr::lead(task) & ideology.int>=2 & dplyr::lead(ideology.int)==1) |
      (task==lag(task) & ideology.int>=3 & lag(ideology.int)==2) |
      (task==dplyr::lead(task) & ideology.int>=3 & dplyr::lead(ideology.int)==2) |
      (task==lag(task) & ideology.int>=4 & lag(ideology.int)==3) |
      (task==dplyr::lead(task) & ideology.int>=4 & dplyr::lead(ideology.int)==3)
  )

## leader with lower political experience (politics > business > media)
# france
dat.e1fr.lessleadexp <- dat.e1fr %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & leadexp=="Business" & lag(leadexp)=="Politics") |
      (task==dplyr::lead(task) & leadexp=="Business" & dplyr::lead(leadexp)=="Politics") |
      (task==lag(task) & leadexp=="Media" & lag(leadexp)=="Politics") |
      (task==dplyr::lead(task) & leadexp=="Media" & dplyr::lead(leadexp)=="Politics") |
      (task==lag(task) & leadexp=="Media" & lag(leadexp)=="Business") |
      (task==dplyr::lead(task) & leadexp=="Media" & dplyr::lead(leadexp)=="Business") 
  )
# Israel
dat.e1il.lessleadexp <- dat.e1il %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & leadexp=="Business" & lag(leadexp)=="Politics") |
      (task==dplyr::lead(task) & leadexp=="Business" & dplyr::lead(leadexp)=="Politics") |
      (task==lag(task) & leadexp=="Media" & lag(leadexp)=="Politics") |
      (task==dplyr::lead(task) & leadexp=="Media" & dplyr::lead(leadexp)=="Politics") |
      (task==lag(task) & leadexp=="Media" & lag(leadexp)=="Business") |
      (task==dplyr::lead(task) & leadexp=="Media" & dplyr::lead(leadexp)=="Business") 
  )
# Sweden
dat.e1sw.lessleadexp <- dat.e1sw %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & leadexp=="Business" & lag(leadexp)=="Politics") |
      (task==dplyr::lead(task) & leadexp=="Business" & dplyr::lead(leadexp)=="Politics") |
      (task==lag(task) & leadexp=="Media" & lag(leadexp)=="Politics") |
      (task==dplyr::lead(task) & leadexp=="Media" & dplyr::lead(leadexp)=="Politics") |
      (task==lag(task) & leadexp=="Media" & lag(leadexp)=="Business") |
      (task==dplyr::lead(task) & leadexp=="Media" & dplyr::lead(leadexp)=="Business") 
  )

## loss of party incumbency (government > outsider > opposition)
# france
dat.e1fr.lessstatus <- dat.e1fr %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & status=="New party" & lag(status)=="Government") |
      (task==dplyr::lead(task) & status=="New party" & dplyr::lead(status)=="Government") |
      (task==lag(task) & status=="Opposition" & lag(status)=="Government") |
      (task==dplyr::lead(task) & status=="Opposition" & dplyr::lead(status)=="Government") |
      (task==lag(task) & status=="Opposition" & lag(status)=="New party") |
      (task==dplyr::lead(task) & status=="Opposition" & dplyr::lead(status)=="New party") 
  )
# Israel
dat.e1il.lessstatus <- dat.e1il %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & status=="New party" & lag(status)=="Government") |
      (task==dplyr::lead(task) & status=="New party" & dplyr::lead(status)=="Government") |
      (task==lag(task) & status=="Opposition" & lag(status)=="Government") |
      (task==dplyr::lead(task) & status=="Opposition" & dplyr::lead(status)=="Government") |
      (task==lag(task) & status=="Opposition" & lag(status)=="New party") |
      (task==dplyr::lead(task) & status=="Opposition" & dplyr::lead(status)=="New party") 
  )
# Sweden
dat.e1sw.lessstatus <- dat.e1sw %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & status=="New party" & lag(status)=="Government") |
      (task==dplyr::lead(task) & status=="New party" & dplyr::lead(status)=="Government") |
      (task==lag(task) & status=="Opposition" & lag(status)=="Government") |
      (task==dplyr::lead(task) & status=="Opposition" & dplyr::lead(status)=="Government") |
      (task==lag(task) & status=="Opposition" & lag(status)=="New party") |
      (task==dplyr::lead(task) & status=="Opposition" & dplyr::lead(status)=="New party") 
  )

## lower female representation 
# france
dat.e1fr.lesswomen <- dat.e1fr %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & women=="Below average" & lag(women)=="Above average") |
      (task==dplyr::lead(task) & women=="Below average" & dplyr::lead(women)=="Above average")
  )
# Israel
dat.e1il.lesswomen <- dat.e1il %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & women=="Below average" & lag(women)=="Above average") |
      (task==dplyr::lead(task) & women=="Below average" & dplyr::lead(women)=="Above average")
  )
# Sweden
dat.e1sw.lesswomen <- dat.e1sw %>% 
  group_by(Response.ID, task) %>%
  filter(
    (task==lag(task) & women=="Below average" & lag(women)=="Above average") |
      (task==dplyr::lead(task) & women=="Below average" & dplyr::lead(women)=="Above average")
  )

### estimate models (marginal means)
# France
out.mm.full.fr1 <- cj(data = dat.e1fr,
                      selected ~ ideology + selection + leadexp + status + women,
                      id = ~Response.ID,
                      estimate = "mm",
)
out.mm.lessideology.fr1 <- cj(data = dat.e1fr.lessideology,
                              selected ~ ideology + selection + leadexp + status + women,
                              id = ~Response.ID,
                              estimate = "mm",
)
out.mm.lessleadexp.fr1 <- cj(data = dat.e1fr.lessleadexp,
                             selected ~ ideology + selection + leadexp + status + women,
                             id = ~Response.ID,
                             estimate = "mm",
)
out.mm.lessstatus.fr1 <- cj(data = dat.e1fr.lessstatus,
                            selected ~ ideology + selection + leadexp + status + women,
                            id = ~Response.ID,
                            estimate = "mm",
)
out.mm.lesswomen.fr1 <- cj(data = dat.e1fr.lesswomen,
                           selected ~ ideology + selection + leadexp + status,
                           id = ~Response.ID,
                           estimate = "mm",
)
# Israel
out.mm.full.il1 <- cj(data = dat.e1il,
                      selected ~ ideology + selection + leadexp + status + women,
                      id = ~Response.ID,
                      estimate = "mm",
)
out.mm.lessideology.il1 <- cj(data = dat.e1il.lessideology,
                              selected ~ ideology + selection + leadexp + status + women,
                              id = ~Response.ID,
                              estimate = "mm",
)
out.mm.lessleadexp.il1 <- cj(data = dat.e1il.lessleadexp,
                             selected ~ ideology + selection + leadexp + status + women,
                             id = ~Response.ID,
                             estimate = "mm",
)
out.mm.lessstatus.il1 <- cj(data = dat.e1il.lessstatus,
                            selected ~ ideology + selection + leadexp + status + women,
                            id = ~Response.ID,
                            estimate = "mm",
)
out.mm.lesswomen.il1 <- cj(data = dat.e1il.lesswomen,
                           selected ~ ideology + selection + leadexp + status,
                           id = ~Response.ID,
                           estimate = "mm",
)
# Sweden
out.mm.full.sw1 <- cj(data = dat.e1sw,
                      selected ~ ideology + selection + leadexp + status + women,
                      id = ~Response.ID,
                      estimate = "mm",
)
out.mm.lessideology.sw1 <- cj(data = dat.e1sw.lessideology,
                              selected ~ ideology + selection + leadexp + status + women,
                              id = ~Response.ID,
                              estimate = "mm",
)
out.mm.lessleadexp.sw1 <- cj(data = dat.e1sw.lessleadexp,
                             selected ~ ideology + selection + leadexp + status + women,
                             id = ~Response.ID,
                             estimate = "mm",
)
out.mm.lessstatus.sw1 <- cj(data = dat.e1sw.lessstatus,
                            selected ~ ideology + selection + leadexp + status + women,
                            id = ~Response.ID,
                            estimate = "mm",
)
out.mm.lesswomen.sw1 <- cj(data = dat.e1sw.lesswomen,
                           selected ~ ideology + selection + leadexp + status,
                           id = ~Response.ID,
                           estimate = "mm",
)

### plot outcomes (Figure 4; showing only MMS for candidate selection in constrained subsets)
## append outcomes of all countries by constrained variable
# generate country indicators
out.mm.full.fr1$cntry <- "France"
out.mm.lessideology.fr1$cntry <- "France"
out.mm.lessleadexp.fr1$cntry <- "France"
out.mm.lessstatus.fr1$cntry <- "France"
out.mm.lesswomen.fr1$cntry <- "France"
out.mm.full.il1$cntry <- "Israel"
out.mm.lessideology.il1$cntry <- "Israel"
out.mm.lessleadexp.il1$cntry <- "Israel"
out.mm.lessstatus.il1$cntry <- "Israel"
out.mm.lesswomen.il1$cntry <- "Israel"
out.mm.full.sw1$cntry <- "Sweden"
out.mm.lessideology.sw1$cntry <- "Sweden"
out.mm.lessleadexp.sw1$cntry <- "Sweden"
out.mm.lessstatus.sw1$cntry <- "Sweden"
out.mm.lesswomen.sw1$cntry <- "Sweden"
# append outcome dataframes
out.mm.full <- rbind(out.mm.full.fr1,out.mm.full.il1,out.mm.full.sw1)
out.mm.lessideology <- rbind(out.mm.lessideology.fr1,out.mm.lessideology.il1,out.mm.lessideology.sw1)
out.mm.lessleadexp <- rbind(out.mm.lessleadexp.fr1,out.mm.lessleadexp.il1,out.mm.lessleadexp.sw1)
out.mm.lessstatus <- rbind(out.mm.lessstatus.fr1,out.mm.lessstatus.il1,out.mm.lessstatus.sw1)
out.mm.lesswomen <- rbind(out.mm.lesswomen.fr1,out.mm.lesswomen.il1,out.mm.lesswomen.sw1)
# change factor order for plot
out.mm.full$cntry <- factor(out.mm.full$cntry, levels=c("Sweden","Israel","France"))
out.mm.lessideology$cntry <- factor(out.mm.lessideology$cntry, levels=c("Sweden","Israel","France"))
out.mm.lessleadexp$cntry <- factor(out.mm.lessleadexp$cntry, levels=c("Sweden","Israel","France"))
out.mm.lessstatus$cntry <- factor(out.mm.lessstatus$cntry, levels=c("Sweden","Israel","France"))
out.mm.lesswomen$cntry <- factor(out.mm.lesswomen$cntry, levels=c("Sweden","Israel","France"))
## generate separate plots per constrained variable
# full sample for baseline comparison
plot.full <- ggplot(filter(out.mm.full,feature=="selection"),aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.75), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Full Sample\n(No Trade-offs)", x = "Marginal Mean") +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
  )
# loss of ideology
plot.lessideology <- ggplot(filter(out.mm.lessideology,feature=="selection"),aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.75), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Ideological\nDistance", x = "Marginal Mean") +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
  )
# loss of leader experience
plot.lessleadexp <- ggplot(filter(out.mm.lessleadexp,feature=="selection"),aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.75), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Leader\nExperience", x = "Marginal Mean") +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
  )
# loss of party status
plot.lessstatus <- ggplot(filter(out.mm.lessstatus,feature=="selection"),aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.75), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Party\nIncumbency", x = "Marginal Mean") +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
  )
# loss of female representation
plot.lesswomen <- ggplot(filter(out.mm.lesswomen,feature=="selection"),aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.75), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Female\nRepresentation", x = "Marginal Mean") +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
  )
## combine plots
p1 <- plot.full + plot.lessideology + plot.lessleadexp
p2 <- plot_spacer() + plot.lessstatus + plot.lesswomen + plot_spacer() + plot_layout(widths = c(1,3,3,1.75))
plot.cost <- p1 / p2 + plot_layout(nrow = 2, guides = "collect") & theme(legend.position = "bottom")
# save file (Figure 4)
ggsave(file = "Figure4.pdf", plot=plot.cost, width=8, height=6.5)


######################################################
# estimation - importance - full results (Figure B3) #
######################################################
# add missing female representation obs to lesswomen model outcomes
women.filler <- out.mm.lesswomen
women.filler.1 <- data.frame(outcome="selected", statistic="mm", feature="women", level="Below average",
                             estimate=NA, std.error=NA, z=NA, p=NA, lower=NA, upper=NA, cntry="France")
women.filler.2 <- data.frame(outcome="selected", statistic="mm", feature="women", level="Above average",
                             estimate=NA, std.error=NA, z=NA, p=NA, lower=NA, upper=NA, cntry="France")
women.filler.3 <- data.frame(outcome="selected", statistic="mm", feature="women", level="Below average",
                             estimate=NA, std.error=NA, z=NA, p=NA, lower=NA, upper=NA, cntry="Israel")
women.filler.4 <- data.frame(outcome="selected", statistic="mm", feature="women", level="Above average",
                             estimate=NA, std.error=NA, z=NA, p=NA, lower=NA, upper=NA, cntry="Israel")
women.filler.5 <- data.frame(outcome="selected", statistic="mm", feature="women", level="Below average",
                             estimate=NA, std.error=NA, z=NA, p=NA, lower=NA, upper=NA, cntry="Sweden")
women.filler.6 <- data.frame(outcome="selected", statistic="mm", feature="women", level="Above average",
                             estimate=NA, std.error=NA, z=NA, p=NA, lower=NA, upper=NA, cntry="Sweden")
out.mm.lesswomen <- rbind(women.filler,women.filler.1,women.filler.2,women.filler.3,women.filler.4,women.filler.5,
                          women.filler.6)
# change factor order for graphs
out.mm.lessideology$cntry <- factor(out.mm.lessideology$cntry, levels=c("Sweden","Israel","France"))
out.mm.lessleadexp$cntry <- factor(out.mm.lessleadexp$cntry, levels=c("Sweden","Israel","France"))
out.mm.lessstatus$cntry <- factor(out.mm.lessstatus$cntry, levels=c("Sweden","Israel","France"))
out.mm.lesswomen$cntry <- factor(out.mm.lesswomen$cntry, levels=c("Sweden","Israel","France"))
# define attribute order for graphs
out.mm.lessideology$feature <- factor(out.mm.lessideology$feature, 
                                      levels = c("selection","ideology","leadexp","status","women"))
out.mm.lessleadexp$feature <- factor(out.mm.lessleadexp$feature, 
                                     levels = c("selection","ideology","leadexp","status","women"))
out.mm.lessstatus$feature <- factor(out.mm.lessstatus$feature, 
                                    levels = c("selection","ideology","leadexp","status","women"))
out.mm.lesswomen$feature <- factor(out.mm.lesswomen$feature, 
                                   levels = c("selection","ideology","leadexp","status","women"))
# define labels for graphs
attributes <- c("Ideological \nDistance","Candidae \nSelection","Leader \nExperience","Party \nStatus",
                "Female \nRepresentation") 
names(attributes) <- c("ideology","selection","leadexp","status","women")
# loss of ideology
plot.lessideology.all <- ggplot(out.mm.lessideology,aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.8), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Ideological\nProximity", x = "Marginal Mean") +
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
             labeller = labeller(.rows = attributes)) +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    title =element_text(size=8),
  )
# loss of leader experience
plot.lessleadexp.all <- ggplot(out.mm.lessleadexp,aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.8), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Leader\nExperience", x = "Marginal Mean") +
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
             labeller = labeller(.rows = attributes)) +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    title =element_text(size=8),
  )
# loss of party status
plot.lessstatus.all <- ggplot(out.mm.lessstatus,aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.8), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Party\nIncumbency", x = "Marginal Mean") +
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
             labeller = labeller(.rows = attributes)) +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    title =element_text(size=8),
  )
# loss of female representation
plot.lesswomen.all <- ggplot(out.mm.lesswomen,aes(x=estimate,y=level,group=cntry, color=cntry)) +
  geom_point(aes(shape=cntry), position = ggstance::position_dodgev(height = 0.5)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper),  
                 size = 0.6, height = 0, na.rm = TRUE,
                 position = ggstance::position_dodgev(height = 0.5)) +
  geom_vline(aes(xintercept=0.5),linetype="dotted") +
  scale_color_manual(values = c("black","gray40","gray65"), breaks = c("France","Israel","Sweden")) +
  scale_shape_manual(values = c(16,15,17), breaks = c("France","Israel","Sweden")) +
  scale_x_continuous(limits = c(0,0.8), breaks = seq(0,0.75,0.25), labels=c("0","0.25","0.5","0.75")) +
  labs(title = "Trade-off With Female\nRepresentation", x = "Marginal Mean") +
  facet_wrap(~ feature, ncol = 1L, scales = "free_y", strip.position = "right",
             labeller = labeller(.rows = attributes)) +
  theme_bw() +
  ggplot2::theme(
    legend.position="bottom",
    legend.title = element_blank(),
    panel.grid = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    title =element_text(size=8),
  )
## combine plots
plot.lessideology.all + plot.lessleadexp.all + plot.lessstatus.all + plot.lesswomen.all + 
  plot_layout(nrow = 1, guides = "collect") & theme(legend.position = "bottom")

plot.cost.full <- plot.lessideology.all + plot.lessleadexp.all + plot.lessstatus.all + plot.lesswomen.all + 
  plot_layout(nrow = 1, guides = "collect") & theme(legend.position = "bottom")

# save file (Figure B3)
ggsave(file = "FigureB3.pdf", plot=plot.cost.full, width=8, height=8)
